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ABSTRACT 

We simulate mergers between galaxies containing collisionally-relaxed nuclei around massive black 
holes (MBHs). Our galaxies contain four mass groups, representative of old stellar populations; a 
primary goal is to understand the distribution of stellar-mass black holes (BHs) after the merger. 
Mergers are followed using direct-summation A'^-body simulations, assuming a mass ratio of 1:3 and 
two different orbits. Evolution of the binary MBH is followed until its separation has shrunk by a 
factor of 20 below the hard-binary separation. During the galaxy merger, large cores are carved out 
in the stellar distribution, with radii several times the influence radius of the massive binary. Much of 
the pre-existing mass segregation is erased during this phase. We follow the evolution of the merged 
galaxies for approximately three, central relaxation times after coalescence of the massive binary; both 
standard, and top-heavy, mass functions are considered. The cores that were formed in the stellar 
distribution persist, and the distribution of the stellar-mass black holes evolves against this essentially 
fixed background. Even after one central relaxation time, these models look very different from the 
relaxed, multi-mass models that are often assumed to describe the distribution of stars and stellar 
remnants near a massive BH. While the stellar BHs do form a cusp on roughly a relaxation time-scale, 
the BH density can be much smaller than in those models. We discuss the implications of our results 
for the EMRI problem and for the existence of Bahcall-Wolf cusps. 

Subject headings: Galaxy:center - stellar dynamics 



1. INTRODUCTION 

The massive black holes (MBHs) that reside at the 
centers of some nearby galaxies are believed to grow to- 
gether with their hosts through mergers: MBHs grow 
partly as a result of gas accretion, and partly by coales- 
cence with other MBHs that arc brought into the nucleus 
during the merger process (Begclman ct al. 1980). The 
detailed assembly history of MBHs is poorly understood; 
major uncertainties include the "seed" mass distribution 
of MBHs at high redshift, the typical gas accretion effi- 
ciency, and the frequency with which MBHs are ejected 
due to gravitational- wave r ecoil (jKauffmann fc Haehnelti 
120001: iVolonteri eralll2003[ ). But a robust prediction of 
the hierarchical models is that galaxies hosting MBHs in 
the nearby universe were formed from less massive sys- 
tems, at least some of which already contained MBHs. 

Binary MBHs created during galaxy mergers leave im- 
prints on the stellar distribution: for instance, they cre- 
ate low-density cores, by excha nging energy with pass- 
ing s tars (jBegelman et al.iriQSOHMilosavlievic fc MerrittI 
|2001| ) . Such cores are observed to b e ubiquitous in stellar 
spheroids brighte r than ^ IO^^Lq (jFerrarese et al.iri994l : 
iLauer et al.lll995l ) and their sizes ~ of order the influence 
radius of the (presumably single) MBH - are consistent 
with the pre dictions of merger models (jGrahaml 120041 : 
iMerrittI [20060 . Here we define the influence radius as 



rh 



GM, 



(1) 



|alessi ag@nipa-garching.mpg.dc 

'"Rlax-Planck institut fiir Astrophysik, Karl-Schwarzschild- 
Str. 1, D-85748 Garching, Germany 

^ Department of Physics and Center for Computational Rel- 
ativity and Gravitation, Rochester Institute of Technology, 54 
Lomb Memorial Drive, Rochester, NY 



where a is the rms velocity of stars in any direction at 
r <; 7'h. Cores of radius ~ rh become difficult to resolve 
in galaxies beyond the Local Group if the MBH mass 
is below ^ 10^ Mq. Even in the nearest nucleus, that of 
the Milky Way, the presence of a parsec-scale core around 
Sgr A* was only clear ly established in the last few years 
(jlJuchholz et aij[20Q9l ). 

In the absence of MBHs, mergers tend to preserve 
the form of the stellar distribution near the centers of 
galaxies (jDehnenI I2005D . Binary MBHs, however, are 
efficient at er asing the structure t hat was present on 
scales < rh (jMerritt fc Cruzl [20011 ). and this fact pre- 
cludes drawing definite conclusions about the nuclear 
properties of the galaxies that preceded the hosts of ob- 
served MBHs. On the other hand, it is well established 
that low-luminosity galaxies ha ve higher central d ensities 
than high- luminosity galaxies (jKormendvl 119851) . This 
is true both in terms of the mean density within the 
half-mass radius, and also in terms of the density on the 
smallest resolvable scales: low-luminosity spheroids often 
contain dense, nuclear stars clusters (NSCs), with sizes 
of order 10 pc and masses ^1% the mass of the galaxy 
(|BoekeHl2010( ). NSC masses are therefore comparable to, 
or som ewhat greater than, MBH masses (jFerrarese et al.l 
I2006al) . although NSCs have been sho wn to coexist with 
MBHs in only a handfu l of galaxies (jSeth et al.l 120081: 
iGraham fc Spitleil [20091 ) . The NSC in the Milky Way 
is believed to be a representative example: its half-light 
radius is 3 — 5 pc, or 1 — 2rh, a nd its mass is a few times 
10^M(D or s everal times M. (jGraham fc Spitleil [20091 : 
ISch6dell[20Tl . 

In the NSC of the Milky Way, the two-body relaxation 
time is comparable with the age of the universe, and 
this is consistent with the persistence of a core (as op- 
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posed to alBa.hcall fc Wolj (|19760 cusp) in the late- type 
stars (|Merrittil2010D . But in fainter systems, central re- 
laxation times are shorter. For instance, in the Virgo 
cluster, galaxies with NSCs have nuclear half-light re- 
laxation times that scale with host-galaxy luminosity as 
(fMerritt 2009i) 
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By comparison, the mean time between "major mergers" 
(mergers with mass ratios 3 : 1 or less) of dark-matter 
haloes in the hierarchical models varies from '^ 0.2 Gyr 
at redshift z = 10 to '^ 10^° yr at z = 1, with a weak de- 
pendence on halo mass (jFakhouri et al.|[2010[) . This com- 
parison suggests that the progenitors of many spheroids 
in the current universe may have been galaxies contain- 
ing nuclei that were able to attain a coUisionally-relaxed 
state before the merger that formed them took place. 

In the absence of a MBH, collisional relaxation implies 
mass segregation and core collapse. If a MBH is present, 
mass segregation still occurs, but core collapse is inhib- 
ited by the fixed potential due to the MBH. A collisional 
steady state, which is reached by a time ~ Tr(rh) at radii 
r ^ rh, is characterized by a Bahcall-Wolf, n ~ r~^/'* 
cusp in the dominant component at r <^ 0.2rh. If there 
is a mass spectrum, less-massive objects follow a shal- 
lower profile, n ~ r~^^^, while mo re-massive objects fol- 
low a steep er profile, n ~ r~^ (jBahcall fc WolJ 119771 : 
iHopman fc Alexanderll2006l ). 

As a first approximation, the mass spectrum of an 
evolved stellar population can be represented in terms 
of just two components: objects of roughly one solar 
mass or less (main-sequence stars, white dwarfs, neu- 
tron stars); and remnant black holes (BHs) with masses 
10 — 20 Mq. Standard initial mass functions predict 
that roughly 1% of the total mass will be in stellar BHs 
(I Alexandeii I2005D : so-ca lled "top-heavy" mass functions 
(e.g. 'Bart ko et ai]|2010D predict a larger fraction. In a 
collisionally relaxed nucleus, the density of stellar BHs 
will rise more steeply toward the center than the density 
of the stars. Mergers between galaxies with such nuclei 
would be expected to modify these steady-state distribu- 
tions substantially, and also to affect (increase) the time 
scale over which a collisionally r elaxed cusp could be re- 
generated following the merger (iMerritt fc Szelll [20061 ). 

These arguments motivated us to carry out merger sim- 
ulations between galaxies containing multi-component, 
mass-segregated nuclei a round MBHs. As in previous 
papers from this s e ries (iMilosavlievic fc MerrittI 120011 : 
iBerczik et all l2005t IMerritt et al.l l2007b[ ). our merger 
simulations are purely stellar-dynamical. In some galax- 
ies, torques from gas would assist in the evolution 
of binary MBHs (lEscala et al.l 120051: iDotti et al.l 120071 : 
iCuadra et al.l 120091 )^ Gas also implies star formation, 
and there is evidence for complex star for mation his- 
tories in many NSCs (jWalcher et al.ll2006[) . But TV- 
body simulations that allo w for non-sphe rical geome- 
tries (Bcrczik et al. 2006; Khan et al.ll201l[ ) have shown 
that purely dissipationless energy exchange with ambient 
stars can bring binary MBHs to milliparsec separations 
on time scales much shorter than galaxy lifetimes. Un- 
less the late evolution of the binary MBH is greatly ac- 
celerated by torques from the gas, the influence of the bi- 



nary on the distribution of the stellar populations should 
be accurately reproduced by our dissipationless models. 
With respect to star formation, population synthesis of 
NSC spectra suggest that most of the mass typically re- 
sides in stars with ages of order 5-10 Gyr (jFiger et al.l 
2004; Boeker 20io[ ). i.e. an old population. 

Simulating the merger of galaxy-sized systems, while 
enforcing the spatial and temporal resolution required to 
faithfully reproduce the dynamics of stars on scales <iC rh 
around the central MBH, is computationally demand- 
ing. Our simulations used '~ 10^ particles per galaxy, 
and the models were adva nced using a paral lel, direct- 
summation iV-body code (iHarfst et al.ll2007D . The in- 
tegrations were accelerated using special-purpose hard- 
ware. The galaxy models contained four mass groups, 
representing an evolved stellar population. Initial condi- 
tions of the merging galaxies were constructed in a two- 
stage process: models of mass-segregated nuclei around 
a MBH were first constructed, then these collisionally- 
relaxed models were imbedded into larger, spheroid- sized 
models. Mergers were then carried out, assuming a 
galaxy mass ratio of 1 : 3. 

A major motivation for our new simulations was 
the need to understand the distribution of stellar rem- 
nants, particularly stellar-mass BHs, near the centers 
of galaxies. Knowledge of the BH density well inside 
r\i is crucial for predicting the rates of many astro- 
physically interesting processes; in particular, the rate 
of capture of stellar-mass BHs by MBHs, or EMRIs 
(jAmaro-Seoane et al1l2007[ ). Published EMRI rate calcu- 
lations almost always assume a state of mass segregation, 
implying a high density of ste l lar remnants ne ar the MBH 
(jHopman fc Alexandeii [2006t lHopma3l2009D . However, 
in a nucleus formed via a merger, any pre-existing mass 
segregation would have been disrupted by the binary 
MBH when it created a core; whether or not the massive 
remnants would have had time to re-segregate following 
the merger is difficult to assess without full iV-body sim- 
ulations. 

Our simulations followed the evolution of the binary 
MBHs for a time almost long enough that gravitational 
wave emission would dominate the binaries' evolution. 
We then combined the two MBH particles into one, sim- 
ulating gravitational wave coalescence, and continued the 
A^-body integrations for a time corresponding to several 
relaxation times at the (new) infiuence radius. In this 
post-merger evolutionary phase, we also considered the 
consequences of varying the relative numbers of the dif- 
ferent mass components. In this way, we were able, for 
the first time, to observe how rapidly the stellar BHs 
would re-segregate following a merger. We found sub- 
stantially longer time scales for this evolution than in 
earlier simulations that started from physically less mo- 
tivated initial conditions. 

The paper is organized as follows. In Section [2] we de- 
scribe the procedure to generate equilibrium segregated 
models starting from single-component models. Scaling 
to physical units is discussed in Section |3| Evolution of 
the binary MBH and its effects on the underlying stel- 
lar distribution are described in Section |4| In Section [5] 
we describe the evolution of the light and heavy objects 
after the massive binary has undergone coalescence. Sec- 
tion |6| describes the shapes and kinematics of the merger 
remnants. Section |7| discusses the implications of our re- 
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suits for the formation and observability of Bahcall-Wolf 
cusps, and for the distribution of stellar remnants. 

2. INITIAL MODELS AND NUMERICAL METHODS 

Multi-mass Fokker-Planck models have been 
constructed for stars around a MBH (e.g. 
iHopman fc Alexandeii |2006[) . Extending these models 
beyond rh - where the stellar distributions are ex- 
pected to be unrelaxed and where the gravitational 
potential contains contributions from stars as well as 
from the MBH - is problematic. Instead, we used 
TV-body integrations to create models of galaxies with 
coUisionally-relaxed nuclei. The major difficulty was 
obtaining sufficient resolution on the scale of the relaxed 
density cusp, without using a prohibitively large number 
of particles overall. In brief, we proceeded as follows. 

1. A mass-segregated model of the inner parts of a 
galaxy containing a MBH was created via A^-body in- 
tegrations, starting from a configuration in which the 
different mass groups all had the same phase-space dis- 
tribution. This model had a total mass of 50M,; scaled 
to a galaxy like the Milky Way, the outer radius of the 
model would be ^ 10 pc. 

2. Smooth representations of the density profiles were 
constructed for each of the A^-body species at the end 
of the integration. These functions were then "spliced" 
onto a larger, unrelaxed model, at a radius where the 
effects of mass segregation were essentially zero. This 
larger model had a mass of 200 Af,, or roughly 1/5 the 
mass of an entire stellar spheroid. 

3. The smooth functions representing the different 
mass components in this larger model were used to gen- 
erate Monte-Carlo positions and velocities as initial con- 
ditions for the A^-body integrations. 

4. Two such A^-body models, with different total 
masses and radii, were placed in orbit around each other 
and integrated forward until the two MBH particles had 
formed a tightly-bound pair at the center. At this time, 
the merged galaxy contained a large, low-density core 
created by the binary MBH. 

5. The two MBH particles were combined into a single 
particle and the merged galaxy was re-sampled using a 
smaller A^. This model was then integrated forward for a 
few central relaxation times, allowing the different mass 
groups to again evolve toward a collisional steady state 
near the center. 

In more detail, the mass-segregated models described 
in step 1 were created as follows. 

Initial conditions were generated from a density law 
having roughly the expected, steady-state distribution 
near the MBH. We used the modified Prugnie l- Simie n 
(1997) model described by iTerzic fc GrahamI (|2005f ). 
which has a central density cusp of adjustable slope: 



pir)=p 
X exp 



l+(^ 
r 



-y/a 



{-M(r"+r-)/rgs]i/""}. (3) 

Setting 7 = 3/2 in the first term gives p ^ r~^/^ near 
the center, which is close to the collisionally-relaxed 
density profile expected for the dominant populatio n 
in a multi-mass cusp fe.g. iHopman fc Alexandeii I2006D . 
The parameter Vg determines the extent of the cusp; 




Fig. 1. — Thin, dotted (black) curve is the density profile of 
the initial, non-mass-segregated model, equation Q. Other curves 
show mass density profiles of the four species after integration for 
~ 1.5 central relaxation times: 1 Mq main-sequence stars (thin, 
red); 0.6 Mq white dwarfs (dashed, green); 1.6 Mq neutron stars 
(dash-dotted, blue) and 10 Mq stellar black holes (thick, black). 
Each curve represents the combined density of eight independent 
integrations with 132k particles. Scaled to the Milky Way, the 
unit of length is approximately 10 pc and the total mass is approx- 
imately 2 X 10* Mq. Thick dotted (black) curve shows the second, 
more extended analytical model into which the mass-segregated 
cusp was imbedded; this model has four times the mass of the first 
model. 

in relaxed, single-com ponent models, r^ w 0.2rh (e.g. 
iMerritt fc Szellll2006D . The parameter p sets the power- 
law slope beyond the central cusp; we set p = 0.5, i.e. a 
relatively constant density like that of the nuclear stellar 
disk of the Alilky Way. The two final terms on the right 
hand side of equation ([3]) mimic a deprojected Sersic-law 
galaxy, with n the Sersic index; the parameter b can be 
related to n if rps is ide ntified with the effective radius 
(iTerzic fc GrahamI |2005( ) . In our case, the exponential 
term is invoked only to provide a sharp truncation to the 
model outside of ~ a few rh; we set n — 1.5. Given these 
parameters, and setting a = 4, we could then solve for 
the cusp radius r^ in units of the model scale- length rpg : 
rg w 0.05rps. The result is shown as the thin dotted 
curve in Figure[T] 

This initial model was assigned a mass of 50M,. If we 
equate A/, and rh with their values in the Milky Way 
(respectively ^ 4 x 10^ Mq and ~ 2.5 pc), the unit of 
length is ^ 10 pc and the total mass is ^ 2 x 10® Mq. 
Henceforth we refer to this subsystem as the NSC. 

The NSC model was assumed to be made up of four 
discrete mass groups. The relative values of the particle 
masses were 1 : 0.6 : 1.4 : 10. These represent, respec- 
tively, one-solar-mass main sequence stars (MS); white 
dwarfs (WD); neutron stars (NS); and 10 M© black holes 
(BH). The relative numbers of the four populations were 
set to 

Nms ■■ NwD ■■ Nns ■■ Nbh = 1 : 0.2 : 0.02 : 0.005, (4) 

independent of radius. In reality, these fractions would 
depend on the initial mass function and the star forma- 
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tion history. Standard values are 

Nms ■■ NwD ■■ Nns ■■ Nbh w 1 : 0.1 : 0.01 : 0.001 (5) 

(e.g. lAlexandeil I2005D . By comparison, our model con- 
tains roughly twice the numbers of WDs and NSs and 
five times the number of BHs. This was done in order 
to improve the statistics for the remnant populations, 
particularly the BHs. Our choices could also be seen as 
corresponding to a. "top-heavy" initial mass function (e.g. 
iManess et al.l 120071 ) . In our post- merger simulations, we 



explored the consequences of varying these ratios. 

Initial positions and velocities of particles from the 
four mass groups were generated in a standard way: (i) 
The gravitational potential $(r) was computed from the 
known mass distribution (equation|3]) plus the central 
MBH particle, (ii) The function 7V(< u, r), the cumu- 
lative distr ibution of velociti es at each radius, was com- 
puted as in lSzell et al.l ([2005) from p{r) and ^{r), assum- 
ing an isotropic velocity distribution, (iii) Monte-Carlo 
positions and velocities were generated from A^(< v, r) 
for each of the four mass groups, with the relative num- 
bers given by equationlH 

We generated eight such models, each containing 
131071 particles, using different initial seeds for the ran- 
dom number generator in each case. We then indepen- 
dently integrated these eight models forward for a time 
corresponding to roughly 1.5 central relaxation times as 
defined by equation ([7]) (using for m the mass of a MS 
particle). The initial relaxation time (which is inde- 
pendent of radius near the MBH in a p ~^ ^-3/2 ^.^gp^ 
was roughly 130 in model units {G — M — rps — 1). 
We used the direc t-summation A^-body code ^GRAPE 
(jHarfst et al.ll2007l ). integrating each model on one node 
of the RIT GRAPE cluster. Figure^] shows the evo- 
lution of the Lagrange radii for the four mass groups. 
This figure combines data from the eight independent 
iV-body integrations; the total number of particles rep- 
resented is roughly one million. The stellar BHs (of 
total number 8 x 635 = 5080) accumulate toward the 
center; by t = 200, their distribution appears to have 
reached a steady state inside ~ 0.1 « Ipc. The lighter 
populations evolve less, as expected, since their distribu- 
tion was close to the steady-state form at the start (by 
design). The final density profiles of the four popula- 
tions are shown in Figure[T] The BHs follow p ^ r~'^ at 
r <, O.l Ri I pc, and they dominate the total (mass) den- 
sity inside r « 0.007 ~ 0.07 pc. The density cusp defined 
by the less massive populations is nearly unchanged aside 
from a slight expansion due to heating by the BHs. 

The mass-segregated model was then imbedded into 
a more extended, unsegregated model with four times 
the total mass. The template for this larger model is 
shown as the thick dotted line in Figure[Il It follows 
the same density law as in equation ([3]) , but with larger 
scale length r'pg. The imbedding was achieved as follows: 
Smoothed density profiles were constructed from each of 
the four mass groups in the evolved A^-body model (the 
same functions plotted in FigureH]). These smooth pro- 
files were then matched onto the density of the extended 
model at a radius ~ 0.3rp5'; beyond this radius, essen- 
tially no evolution occurred in the A^-body integrations 
and so the A'^-body density profiles (with appropriate 
vertical scalings) matched well onto the analytic profile. 




20 40 50 80 100 120 140 160 180 200 
time 

Fig. 2. — Lagrange radii of the four stellar species during the TV- 
body integrations of the initial model shown in Figurc[T] This plot 
is based on the roughly one million particles in the 8, independent 
integrations of that model generated with different random num- 
ber seeds. Red: 1 Mq main-sequence stars; green: 0.6 M© white 
dwarfs; blue: 1.6 Mq neutron stars; black: 10 Mq stellar BHs. The 
BHs form a dense cluster around the MBH in approximately one, 
central relaxation time as defined by the main-sequence stars. The 
lighter populations are pushed slightly outward during this time. 
Scaled to the Milky Way, the approximate unit of length is 10 pc. 

A small degree of smoothing was nevertheless necessary 
near the matching radius to keep the first derivatives of 
the density continuous. Finally, small (a few percent) ad- 
justments were made in the vertical normalizations of the 
four density profiles in order to recover precisely the orig- 
inal ratios between the total numbers of the four species. 

The total mass of this extended model was 200 times 
the MBH mass. Since observed bulge masses are ~ 
10^ Af,, such a model can be interpreted as comprising 
the innermost ~ 20% of a real bulge. In what follows, 
we refer to this model as "the bulge." 

A^-body realizations of the bulge models were then con- 
structed in the same way as described above. Two such 
models were required for each merger simulation. We 
considered unequal-mass mergers with a mass ratio of 
3:1. The radius of the smaller bulge was scaled as the 
square-root of the bulge mass. Particle masses (aside 
from the MBH particle) were the same in the two bulges, 
i.e., particle number scaled linearly with total mass. 

Table[l] lists the important parameters of the merging 
systems: the mass ratio of the two MBHs (equal to the 
galaxy mass ratio), the ratio of MBH mass to host bulge 
mass, the ratio of the MS stars mass to MBH mass, the 
number of particles in each galaxy, and the initial sepa- 
ration between the bulge centers, Ar, expressed in units 
of the outer radius of the larger bulge, i?i. The table 
also gives the influence radius associated with the MBH 
in the larger galaxy, under two definitions. In addition to 
the first definition given in equation[Tl we also compute a 
second, mass-based infiuence radius: the radius r,ii that 
contains a mass in stars equal to twice the MBH mass: 

A'f^r <r„,)-2Af.. (6) 

The mass-based definition is the most straightforward to 
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TABLE 1 
Galaxy merger parameters 



Ml : Ma M, I Afgai Mms/M. Ni N2 Ar/Ri 



rh 



3:1 



0.005 



0.0005 



400k 120k 



1.5 



0.10 



TABLE 2 
Post-merger models 



Model 


N 


Fractions 


M./M,,i 


Al 


250k 


1:0.2:0.02:0.005 


0.0050 


A2 


225k 


1:0.1:0.01:0.001 


0.0055 


Bl 


250k 


1:0.2:0.02:0.005 


0.0050 


B2 


225k 


1:0.1:0.01:0.001 


0.0055 



apply in iV-body models like ours; in computing rh, a 
choice must be made about the radius at which to evalu- 
ate a, and this, combined with the need to bin particles, 
can lead to factor of ~ two uncertainties in r-h- The val- 
ues of Th and Tm in Table[T] refer to the larger of the two 
galaxies. 

The two bulge models were placed far enough apart 
initially that there was no overlap, but not much farther, 
in order to minimize the integration time. Two relative 
orbits were considered: a circular orbit (Model A), and 
an eccentric orbit in which the initial relative velocity 
was set to 0.7 times the circular value (Model B). 

The merger simulations were also carried out using 
0GRAPE, both in combination with GRAPE hard- 
ware at the Rochester Institute of Technology and with 
GPU hardware at the Max-Planck Institute for Astro- 
physics in Garching by means of the Sapporo library 
(jGaburov et al.ll2009f). We conser vatively set the accu- 
racy parameter (jHarfst et al.ll2007l ) to 0.005 and the soft- 
ening to lO"** (in iV-body units). Such small values of 
the softening and accuracy parameters were necessary to 
accurately follow the dynamics of the massive binary. 

The simulations were continued until the two MBH 
particles had formed a tight binary and the binary sep- 
aration had shrunk by an additional factor of ^ 20. At 
this time, the two MBH particles were combined into a 
single particle, which was placed at the center-of-mass 
position and velocity of the binary. A random subset 
of particles was then chosen from this model, yielding 
a new model with the same phase-space distributions of 
the four components but a smaller total N. This "merged 
galaxy" model was then integrated forward, for a time 
corresponding to a few central relaxation times. The sub- 
sampling was a necessary compromise to keep the phys- 
ical integration time from becoming prohibitively long 
while still allowing the simulations to proceed for more 
than one relaxation time. For each of the models pre- 
sented in Table[l] we generated two different submodels 
with different number fractions, as illustrated in Tabled 
In order to generate the models with the standard frac- 
tions Nms ■■ NwD ■■ Nns ■■ Nbh « 1 : 0.1 : 0.01 : 0.001, 
we deleted an appropriate number of particles in each 
mass group from the models with Nms ■ Nwd '■ N^s '■ 
Nbh = 1 : 0.2 : 0.02 : 0.005. 

3. SCALING 



In any A^-body simulation, an important consideration 

is how to relate the computational units of mass, length 

and time to physical units. In our simulations, the mass 

-scale is most naturally set by the mass of the particle (s) 



Q 2tf^p resenting the MBH(s), and the length scale by the 
— TTcdius that encloses a mass in stars that is some multiple 
of the MBH mass. A natural choice for the latter is r^, as 
defined in equation ([6]) . Identifying r^ in the simulations 
with Tm in a real galaxy is only justified, of course, to the 
extent that the radial dependence of the density near the 
MBH is similar in both systems. 

Scaling of the time is more subtle. Two basic time 
scales are of interest: the crossing time, which is deter- 
mined by the total mass and size of the model; and the 
relaxation time, which depends as well on the masses of 
the particles, or equivalently on N: 



Tr 



0.33 cr3 
G^ nw?\n.K 



(7) 



()SpitzeHll987D . Here, n is the stellar number density, m is 
the mass of one star, and In A is the Coulomb logarithm. 
The particle masses in our simulations have the correct 
ratios with respect to one another, but their masses in 
relation to the total galaxy mass is much larger than in 
a real galaxy - due of course to the fact that our N is 
much smaller than 10^^. While there is a well-defined 
relaxation time in the models, that time is much shorter 
relative to the crossing time than it would be in real 
galaxies. 

These statements apply to many iV-body simulations 
that extend over relaxation time scales. A novel feature 
in our simulations is the treatment of the galaxy merger. 
Such a merger requires of order a few crossing times in 
order to reach a (collisionless) steady state. Since cross- 
ing times are much shorter than relaxation times, both 
in reality and in our models, a galaxy would not undergo 
a significant amount of collisional relaxation during the 
merger. During this phase of the simulations, therefore, 
the appropriate unit of time is the crossing time. In the 
phases preceding and following the merger, the appropri- 
ate unit of time is the relaxation time. 

Subtleties arise when one considers the massive bi- 
nary. After the galaxy merger is essentially complete, 
the binary MBH continues to evolve. In a spherical nu- 
cleus, the binary quickly ejects stars on intersecting or- 
bits, and continued hardening takes place on a relax- 
ation time scale, as stars are scattered by other stars 
onto previously-dep leted orbits that intersect the binary 
(jMakino fc Funatoll2004 iMerritt et"aIll2007bD . In such 
models, there is effectively just one time scale - the relax- 
ation time - that determines both the rate of collisional 
evolution of the galaxy and of the central binary follow- 
ing the merger. 

Another mode of evolution is possib le, if the galaxy 
potential is significantly nonspherical (jMerritt fc PoonI 
120041 : iBerczik et al.l I2OO60 In this case, orbital angular 
momenta of stars near the massive binary evolve due 
both to encounters, and to torques from the overall stellar 
potential. Typically the latter dominates, and the sup- 
ply of stars to the massive binary remains high in spite of 
ongoing, slingshot ejections. The binary evolves at a rate 
that is fixed essentially by stellar orbital periods, i.e. by 
the crossing time, and not by the relaxation time. Such 
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evolution appears to be the norm when the galaxy host- 
ing the ma ssive binary was formed in a realisti c merger 
simulation (|Khan et al.ll2011l: iPreto et al.ll2011[ ). and we 
find the same result in our new simulations. This result 
simplifies the time scaling of our models, since it means 
that a single physical time - the crossing time - sets the 
rate of evolution, both during the galaxy merger, and 
immediately afterwards, as the massive binary hardens. 
One potential difficulty does arise, however. If hard- 
ening of the massive binary is simulated for a time that 
is comparable with the A'^-body central relaxation time, 
some collisional evolution in the stellar distribution will 
occur. This may or may not be realistic, since in a real 
galaxy, the ratio of relaxation time to crossing time is 
much larger than in the simulations. For this reason, the 
models that we adopt at the start of the final phase of our 
simulations - single galaxies containing merged MBHs - 
may exhibit overly-segregated nuclei, causing the subse- 
quent, collisional relaxation to occur more quickly than 
it would in a real galaxy. 

4. THE GALAXY MERGER 
4.1. Evolution of the massive binary 

Figures[3] and S] illustrate the large-scale evolution of 
the bulge models described in Table[T] during the galaxy 
merger phase. The two bulge models start out either on 
circular (Model A) or eccentric (Model B) orbits about 
their common center of mass, with initial separations 
roughly one-half the radius of the larger bulge. Due to 
its more eccentric initial orbit. Model B evolves more 
quickly. The trajectories of the MBHs reflect the initial 
orbits of the parent systems, as can be seen in Figure[5j 

Evolution of a binary MBH can be divided roughly into 
four phases. (1) Formation of a bound pair. Merger of 
two galaxies with central MBHs brings the two MBHs 
together, in a time comparable with the galaxy merger 
time, i.e. a few galaxy crossing times. (2) Formation of 
a hard binary. The separation between the two MBHs 
decreases very rapidly, due at first to dynamical friction 
against the stars, and then to the gravitational slingshot: 
near stars are ejected after gaining energy from the mas- 
sive binary. A core is formed at this stage, with size 
roughly equal to the initial separation of the bound pair. 
(3) Binary hardening. If orbits of stars ejected by the 
gravitational slingshot are replenished, the massive bi- 
nary will continue to shrink. Hyper-velocity stars can be 
produced in this process. The stellar core will continue 
to grow. (4) Coalescence. If replenishment of stellar or- 
bits continues, the binary separation decreases to a value 
such that emission of gravitational waves dominates its 
evolution, and the two MBHs coalesce. 

Dynamical friction drives the evolution down to a sepa- 
ration a/ at which the stellar mass enclosed in the binary 
is of order twice the mass of the smaller MBH: 



M{<a-f) «2M2. 



(8) 
M2/M1 than the 



This separation is smaller by a factor ^ 
radius of infiuence of the larger MBH. 

At separations smaller than Of , the binary hardens due 
to gravitational slingshot interactions with passing stars. 
A binary is defined as "hard" when it reaches a separa- 
tion 

an « -j^ (9) 



called the hard-binary separation; a binary is "hard" 
when its binding energy per unit mass, \E\/{AIi + M2), 
exceeds a^. 

The binary enters the gravitational wave (GW) regime 
when the time scale for coalescence due to emission of 
gravity waves: 



T, 



GW- 



256F(e) G3 ^ (Mi + M^f 

5.8xl0"yr/ a \ ViO^Mq \ / IO^Mq " ^ 



F{e) 



0.1 pc 



/i J \Mi + M2 



(10) 



becomes shorter than the time for hardening due to stel- 
lar interactions. Here 



F{e) = (1 - e^ 



-7/2 



73 2 37 

^+24^+96^ 



and n = Ml M2/(Mi -I- M2) is the reduced mass. 

The evolution of the separation between the MBHs is 
shown in FigurelHl The first part of this evolution, which 
is driven by dynamical friction, ends when the separa- 
tion reaches ~ a/, equation ([8]); ay « 0.9 in model units. 
The time to reach this separation was tf « 440 and 
~ 160, respectively, for models A and B. At about the 
same time, gravitational slingshot encounters with the 
stars begin to dominate the binaries' evolution and the 
effect on the galaxy structure starts to become apparent. 
Figure[7| clearly shows a decrease in the central densities 
of both the main-sequence stars and the stellar BHs at 
a time t Ki tf. Similar expansions are observed in the 
white dwarf and neutron star distributions. However, 
the stellar BHs are most affected, due to their higher, 
initial central concentration. 

The orbital semi-major axis and eccentricity of the 
massive binary are shown as functions of time in Fig- 
ure[8l Even in the circular-orbit merger (Model A), the 
massive binary forms with a slightly nonzero eccentricity. 
Subsequent evolution of a and e is driven by interactions 
with stars on intersecting orbits. 

We computed th e time-dependent binary hardening 
rate (|Quinlanlll996[) 



d 
di 



(11) 



by fitting a straight line to a~^{t) in small time intervals 
between the time of binary formation and the end of the 
integration. The results are shown in the bottom panels 
of Figure[9] for both models. The top panels show the 
evolution of 1/a. The hardening rate appears roughly 
constant in time, which suggests that the binaries enter 
the hard phase rather quickly, and then continue harden- 
ing at an approximately constant rate. This behavior is 
a defining property of hard binaries, if the stellar back- 
ground is unchanging, i. e. una ffe cted by th e binary. 

Mikkola & Valtone3 (|1992[ ). iQuinlanI (J1996) and 
ISesana et al., (,2006, ) performed three-body scattering ex- 
periments of circular binaries of varying mass ratio and 
hardness to derive estimates of the hardening rate. They 
provided fitting formulae for the dimensionless parame- 
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Fig. 3. — Snapshots from the early evolution of Model A. The different colors refer to the different mass groups: main sequence (blue), 
white dwarfs (red), neutron stars (green), black holes (black). The MBHs are indicated by full circles. 
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Fig. 4. — Snapshots from the early evolution of Model B. Colors are as in FigurefS] 
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Fig. 5.— Trajectories of the MBHs in Model A (left) and B 
(right) during the galaxy merger phase. Dashed line: trajectory 
of the smaller hole. Solid line; trajectory of the larger hole. The 
initial orbit of the galaxy pair lies in the 2 = plane. 
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Fig. 6. — Separation between the MBH particles as a function of 
time. The horizontal lines indicate af, the approximate separation 
at which dynamical friction ceases to be efficient (equation|8ll , and 
a^, the hard-binary separation (cquation[9]| ■ 

ter H which is related to the hardening rate via 

where p is the stellar mass density and a is the one- 
dimensional velocity dispersion, both assumed constant 
in space and time and unaffected by the presence of the 
binary. The H parameter obtained from these scattering 
experiments is a function of binary mass ratio and hard- 
ness; the latter is defined in terms of Vbin/c, the ratio of 
binary circular speed to field-star velocity dispersion. 

We wish to compare the A^-body results for s with 
the predictions from the scattering experiments. Good 
agreement would imply a "full loss cone," i.e. that the 
rate of interaction of the binary with stars is unaffected 
by slingshot ejections. In spherical galaxy models, binary 
hardening rates fall much below the predictions from 
the scattering experiments, since orbital repopulation is 
driven by t wo-body sc attering, which is a slow process 
for large TV (jMakino fc Fu nato 2 004; .Berczik et al.,,2005; 
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Fig. 7. — Evolution of the Lagrange radii of the main sequence 
stars (left) and stellar BHs (right) during the galaxy merger phase, 
for models A (top) and B (bottom). For clarity, only stars initially 
belonging to the larger galaxy are shown. Formation of the binary 
MBH is reflected in the sudden expansion of the central regions, 
at t Si 500 (Model A) and t ^ 200 (Model B). 
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Fig. 8. — Evolution of the semi-major axis and eccentricity of the 
MBH binary in Model A and B, starting from the time when the 
MBHs are formally bound. 



iMerritt et aIll2007bD . 

We defined the theoretically-expected hardening rate 
to be 

su{a) = H{a) (^) (13) 



■^a J D 

where H is taken from the published scattering experi- 
ments, and p and a are measured in our A'^-body models, 
at a distance D from the MBH. (G = 1 in JV -body units.) 
For J? (a) we adopted the fitting formula of iSesana et alJ 
(|2006[) : 

H = A{\ + alaQ)\ aq = 1.05^^ (14) 

with parameters for a 3 : 1 mass ratio: A — 15.82, 7 = 
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derivable from sca ttering experiments; for instance, 
iSesana et all (|2006[ ) give 



K ^ A{1 + a/ao)'' + B , 



(16) 



Fig. 9. — Evolution of the inverse binary semi-major axis (upper) and the binary hardening rate (lower). The points represent the A'^-body 
results while the dashed lines represent the semi-analytic estimates for the different mass groups, assuming a full loss cone regime and 
evaluating the density and velocity dispersion at D = 0.5 from the binary center of mass. 

—0.95. Because p and a vary with position (and time) 
in the iV-body models, our computed values of s will 
depend on D. We found this dependence to be weak, as 
long as D is not too different (i.e. not more than a factor 
of two) from r^^: both p and a are weakly dependent on 
radius for r « r^. In Figure[9]we plot values of sd for 
D = 0.5, slightly larger than estimated influence radii 
in both models. We find that the iV-body hardening 
rates, s, are quite consistent with the analytic predictions 
S£). The contributions to the binary hardening from the 
different mass groups (assumed to scale in proportion to 
their mass densities) are shown in the bottom panels of 
FigureO The MS stars appear to be responsible for most 
of the binary hardening, followed by the white dwarfs and 
the stellar black holes. 

These results suggest that the massive binaries in our 
simulations are roughly in the "full-loss-cone" regime: 
they harden at a rate that is consistent with the expected 
hardening rate of a binary in an undeple ted field of 
stars. This is in agreement wit h the results of lKhan et al.l 
(|2011h and [Preto^TaD (120 lit ), who also found that the 
stalling that occurs in spherical models is absent in sim- 
ulations that start from a stage preceding merger of the 
two galaxies. Apparently, the non-spherical shapes of 
the merger remnants result in a large population of stars 
on "centrophi lic" orbits: saucer orbits in the axisymmet- 
ric geometry (jSridhar fc Touma"1999') , pyramid orbits in 
the triaxial geometry ([Mcrritt & Vasilicv 2010), etc. We 
discuss the shapes of our models in more detail in Sec- 
tion [6] 

Figured] shows that the eccentricity of the massive bi- 
nary remains roughly constant in both models. In a 
nonrotating galaxy, binary eccentricity is expected to in- 
crease gradually with time: 



de d\n{l/a) 

—r = tC ; 

dt dt 



(15) 



and their Table[3]gives values of A and i? for a 1 : 3 binary 
mass ratio. Stars that encounter the binary in a prograde 
(co-rotating) sense tend to circularize it; Sesana ct al] 
(|2011t ) found that when the fraction of corotating stars 
exceeded ~ 0.7 (as opposed to 0.5 for a nonrotating 
galaxy), binaries tend to circularize. We found that our 
models have roughly this fraction (~ 0.7) of corotating 
stars, consistent with the mild eccentricity growth that 
we see. 

In our simulations, binary hardening rates s are essen- 
tially constant with time. If we assume that this remains 
true beyond the end of our simulations, we can use equa- 
tions (ITU) and (|15p to extrapolate the binary elements 
(a, e) to arbitrarily later times. At some point, gravi- 
tational wave emission will dominate the evolution; cal- 
culating when this happens requires assigning physical 
units to the models. We considered two representative 
scalings, based on assumed MBH masses of 4.0 x 10^ Mq 
and 10^ M0 respectively. Scaling factors for length were 
set at 10 pc and 40 pc respectively, yielding physical val- 
ues of the influence radii of ~ 3pc and ~ 12 pc. As 
discussed in Section [3l the unit of time is determined 
in this (full loss cone) regime by orbital periods, not re- 
laxation times; hence the scaling factor for time, [T], 
is related to the scaling factors for mass and length by 
[T] = ^[LY/{G[M]) and the binary hardening rate in 

physical units is [L\~ \T]~ times the A^-body harden- 
ing rate. 

Given a value for s, and scale factors for mass and 
length, the evolution of the binary semi-major axis at 
late times is determined by 



where K > is a second dimensionless rate, also 



da 



da 
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da 



GW 



9/ N da 
-sa\t) + - 



(17) 
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Fig. 10. — Evolution of the binary elements in the A'^-body phase (symbols) and in the semi-analytic phase (lines). For each model, two 
physical scalings are adopted, as discussed in the text. 



TABLE 3 

Hardening phase parameters 



Model 


Ml 


Tnb 


Thd 


a/ 


e/ 


eg 




(M0) 


(yr) 


(yr) 


(mpc) 






A 


4 X 10** 


2.3 X 10^ 


7.3 X lO'^ 


5.6 


0.23 


8 X 10"'^ 




1 X 10** 


2.7 X lO'^ 


3.7 X lO'^ 


37.2 


0.37 


8 X 10-5 


B 


4 X 10^* 


1.0 X lO'^ 


4.8 X lO'^ 


8.0 


0.58 


4 X 10-5 




1 X 10** 


1.1 X lO'^ 


2.6 X lO'^ 


50.2 


0.60 


2 X 10-* 



where the first term on the right hand side represents 
hardening due to interactions with stars, and the sec- 
ond term represents energy lost to gravitational waves. 
The rate of the latter process depends on e as well as a. 
In extrapolating e beyond the end of the A^-body inte- 
grations, we ignored the effect of the galaxies' rotation 
on the binaries' eccentricity growth and simply applied 
equations (fT5|) and (fTBj) . The rate of change of the orbital 
elements due to GW emission is (|Peterslll964[ ) 



da 



GW 



de 
di 



GW 



64 J^(e) 
5' a^ 
304 eG(e) 



where F{e) is given in equation (jlip . 



G(e) - (1 - e^ 



-5/2 



121 
304* 



(18a) 
(18b) 

(19) 



and 



P 



G= 



MiM2{Mi +M2) 



We numerically solved the coupled equations for the 
evolution of the semi-major axis and eccentricity, start- 



ing from values at the end of the iV-body phase or some- 
what earlier, for each of the two sets of scaling factors. 
The results are shown in FigurefTOl and Table[3l In this 
table, the time Tnb represents the time spent in the 
A^-body hardening phase while Thd represents the time 
from the end of the iV-body integration until coalescence; 
the latter time includes both a stellar-interaction driven, 
and a GW dominated, regime. The total evolution time 
from the beginning of the galaxy merger to full coales- 
cence is given by the sum of the two times. For a Milky 
Way type galaxy undergoing a 3 : 1 merger, this time 
is of the order of 10^ yr for an initial circular orbit and 
6 X 10^ yr for an initially moderately eccentric orbit. The 
table also lists the values of the semi-major axis af and 
eccentricity e/ at the transition between the simulated 
A^-body phase and the semi-analytical phase, as well as 
the value of the eccentricity when the separation reaches 
lOrg = lOGMi/c^. At smaller separations, equations 
like (fT8l) are not valid. 



4.2. Core formation 

The pre-merger galaxies had central density profiles 
that were well approximated as power laws with respect 
to radius. Following the merger, the density on scales 
r < Th is strongly modified (lowered) by the action of 
the massive binary. The process of "core scouring" is 
illustrated in Figure[TTJ which compares the spatial den- 
sity profiles of the different species in the larger galaxy 
at four different times. For this plot, only particles that 
were originally associated with the larger galaxy were 
used. By the time the b inary has become hard, a mass 
deficit (jMilosavlievic et al. 2002) is clearly visible in all 
components on scales significantly larger than Vm- Model 
A shows substantially more evolution of the central den- 
sity at early times; at late times the cores in the two 
models are more similar. 
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Fig. 11. — Spatial density profiles of the different mass groups in 
the larger galaxy. For each species, different lines are for different 
times during the merger phase: at the start of the simulation (solid 
lines), at the time when a ~ af (dashed lines), at the time when 
a ~ a^ (dotted lines) and at a late time when a ^ O.la^ (dashed- 
dotted lines). Vertical lines indicate rm, defined as described in 
the text. 

We estimated core sizes in two ways. First, we fit core- 
Sersic profiles to our models at the end of the merger 
phase, and adopted the resulting values of the break ra- 
dius Rb as estimates of the core radius. The projected 
density profiles of galaxies are globally well fit by the 
Sersic (1968) law: 



/(i?)=/(0)exp{-6(i?/i?e)'/"} 



(20) 



where /(O) is the central intensity, Re is the effective 
half-light radius, and n is a parameter controlling the 
curvature in a log-log plot. The ter m b is not a pa- 
rame ter but a function of n (see e.g. iTerzic fc Grahaiiil 
l2005f ). Deviations from this law appear close to the cen- 
ter, where bright galaxies typically show central light 
deficits with respect to the inward extrapolation of the 



TABLE 4 
Core radii 



Model 


fm 


Rb (MS) 


vc (MS) 


Rt (BH) 


re (BH) 


A 
B 


0.38 
0.36 


2.4 
0.9 


0.8 
0.6 


5.1 
3.5 


0.3 
0.2 



Sersi c law while faint galaxies show central light excesses 
(e.g. iGraham fc Guzma^ [200l iFerrarese etall l2006bl: 
ICote et all 20071). Add ing an inner power-law with slope 
7 (Grahan i'etal]|2003 ) yields 



I{R) 
with 



I'll 



(i?b/i?)"]^/"exp{-6(i?" 



(21) 



/' = lb 2-''/"exp [6 2i/("") {Rb/Re)^'"' 



(22) 



the so-called core-Sersic law. Here, a is an additional pa- 
rameter regulating the transition from the inner power- 
law to the outer Sersic law and i?{,, the break radius, 
marks the distance where the profile changes from one 
regime to the other. 

We also considered a second definition of the core ra- 
dius (HCing 1962), as the projected radius Vc at which the 
surface density falls to one-half its central value. Here, 
"central" was taken to be the value at a projected radius 
of O.lr^. 

Both estimates, computed at the time when a '^ 0.1 au, 
are listed in Table|31 separately for the MS stars and the 
stellar-mass BHs. These radii were computed using all 
the particles from the respective mass groups, without 
regard to the galaxy with which they were originally as- 
sociated. We also give r-^, computed using equation ([6]) ; 
we set Af, = Mi + M2 in that equation, and combined 
together all the stellar species from both galaxies when 
computing A/^. 

In the case of the dominant (MS) mass component, 
Table[4] shows that r^ ~ 2r^ for both Models A and B. 
However, Rb is somewhat larger than Tc- The reason is 
that, due to the flatness of the profile, the radius at which 
the inner profile transitions to the outer profile (i.e. the 
break radius) is much larger than the radius at which 
the density falls to half its central value (i.e. the core 
radius). This is illustrated in Figure[T2] for Model B. 

In what follows, we will adopt re as the more robust 
measure of the core radius. We emphasize that the size of 
the core in the MS stars is larger than the infiuence radius 
of the massive binary in both of our models, whichever 
definition of "core radius" or "influence radius" is used. 
One important consequence, discussed in more detail be- 
low, is that regrowth of a Bahcall-Wolf cusp at radii 
^ Tni following coalescence of the two MBHs leaves the 
core structure essentially unchanged. 

4.3. Stellar ejections 

Stars can become unbound during a galaxy merger 
due to both tidal stripping in the early phases of the 
merger and gravitational slingshot interactions with the 
binary MBHs when this is hard. Unbound stars at the 
end of the merger phase are a result of both mecha- 
nisms. We determined the number of unbound stars in 
the two merger simulations by selecting stars with veloc- 
ity in excess of the local escape speed from the system. 
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TABLE 5 
Fraction of unbound stars 



Model 


/c 


/cl 


/c2 


/cMS 


/cWD 


/cNS 


/cBH 


J^/cTs/A^^Sal 


AfeGs/M. 


A 
B 


0.021 
0.014 


0.005 
0.001 


0.016 
0.013 


0.021 
0.014 


0.020 
0.012 


0.024 
0.016 


0.024 
0.018 


0.018 
0.012 


0.55 
0.45 




Fig. 12. — Projected density profile for MS stars in Model B, 
at the time when a ~ O.la^ (points), with superimposed the best 
fitting core-Sersic model (solid line). The vertical dotted line indi- 
cates Tc while the vertical dashed line indicates R\y. The horizontal 
lines indicate the value of the central density Eo and half the cen- 
tral density. 

At any given time, the escape velocity at a distance r 
from the binary center of mass is Vi^s,tJyr) = -^Z— 2$(r), 
where the gravitational po tential $ can be expressed as 
(jBinnev fc Tremainelll987l equation 2-22) 

$(r) == $(/ < r) + $(r' > r) (23) 

'GM{< r) 



with 

$(r' < r) 

r .In r 

(24) 

having defined M{< r) the mass enclosed within a sphere 
of radius r, and 



-AttG- / pir'y^dr' 
f Jo 



$(r' > r) = 



) = -4ttG / 

Jr 



N 



p{r')r'dr' ^-GY,—. for r, > 
»=i ^' 

(25) 
Table[S]reports the total fraction of unbound stars /e at 
the time when a ~ O.la/i, the fraction of unbound stars 
originally belonging to the first and second galaxy /ei, 
/e2, and the fraction of unbound MS, WD, NS and BHs 
Ims, fwD, Ins, /bh, normalized to the total number of 
stars in their mass group. We find that about 2% of all 
stars are ejected by the end of the merger phase. Model 
A, which has a more gradual evolution, produces more 
unbound stars than Model B. Stars from the different 
mass groups show roughly equal probabilities of being 
ejected. 



We distinguished stars ejected via tidal stripping (TS) 
versus gravitational slingshot (GS) based on the time 
they become unbound: if a star becomes unbound be- 
fore the binary has become hard we consider it ejected 
by TS whereas if it becomes unbound after the binary 
has become hard we consider it ejected by GS. Based on 
this selection, we find that TS is responsible for the ejec- 
tion of about 90% of the escapers. The smaller galaxy 
is the most susceptible to TS, as shown by the fact that 
/e2 ^ /el- Table[5]also lists the mass in stars unbound 
by TS in units of the total stellar mass Afgai and by 
the mass in stars unbound by GS in units of the bi- 
nary mass. High velocity escapers, with velocities as 
high as several times the central escape speed, are mainly 
produced by GS. These stars are analogs of the hyper- 
velocity stars detected in the halo of the Milky Wa 
(e.g. iB rown ct al. 2005, 2006; Gualandris et al. 2005| 
Gualandris fc Portegics Zwart. .2007. : .Gvaramadze et al 
20091 ). Stars unbound bv TS tend to have lower veloci- 
ties, and would be less numerous in a real galaxy, due to 
the deeper potential well. 

5. POST-MERGER EVOLUTION 

After the two MBH particles were combined into one, 
we continued the integrations of the iV-body models 
for several relaxation times. The relaxation time of a 
multi-component system is ill-defined. We are primar- 
ily interested in the time scale for collisional evolution 
of the dominant population, the MS stars. As a sim- 
ple estimate of the relaxation time, we used equation ([7]), 
setting m to the mass mMS of a MS particle, and re- 
placing the number density n by a simple summation, 
nt = riMS + nwD + ^ns + nsH- Other reasonable def- 
initions of Tj. were found to give nearly identical nu- 
merical values, a consequence of the fact that the MS 
stars dominate the total numbers at the radii of inter- 
est, and the fact that the masses of the three major 
groups (MS, NS, WD) are very similar. The result- 
ing estimates of T^, computed at the MBH infiuence 
radius r^ at the beginning of the post-merger phase, 
are given in Table|6l For the Coulomb logarithm we 
used InA = \n{rhcr^ /2GmMs) — ln(M,/2mMs), with 
rfi = GM,/a^ and M, the mass of the merged MBHs. 

Ignoring the infiuence of the other components, a cusp 
in the MS stars is expected to re-form in a time of roughly 
Ti{r^), at radii r <, 0.2rni, after being des troyed by the 
massive binary (e.s. IMerritt fc Szelll 120061 ) . In the case 
of the BHs, their central density should increase more 
rapidly, by a factor ~ rnBH /itt-ms = 10, as they segre- 
gate spatially with respect to the lighter components. If 
the density in the heavier (BH) component should ever 
approach locally the density in the lighter components, 
heating of the light particles by the heavy particles will 
occur, causing the density of the former to decrease. 

Figure[T2] shows the evolution of the enclosed mass in 
each species in the post-merger integrations. Within the 
infiuence sphere, the general trend is for the mass in the 
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Fig. 13. — Time evolution of the stellar mass enclosed in the MBH influence sphere, for each species, in the post-merger models. 



TABLE 6 
Central properties at the start of the post-merger phase 



Model 


^h 


r-m 


In A 


T,{rm) 


Al 


0.09 


0.40 


6.4 


4.7 X 10^ 


A2 


0.09 


0.43 


6.4 


5.8 X 10« 


Bl 


0.09 


0.35 


6.4 


3.8 X 103 


B2 


0.09 


0.37 


6.4 


3.7 X 10» 



lighter coraponents to decrease with time, while the mass 
in BHs increases. The mass in NSs decreases in models 
Al and Bl while it appears nearly constant in models A2 
and B2. In terms of mass density, Figure[T3] shows that 
the BHs quickly (on a time scale <^ T^irm)) form the 
expected, steep density profile, p '^ r^^. However, the 
MS stars maintain a core-like profile. Only at very small 
radii, r < 0.2r^, does evolution toward a cusp appear to 
occur in the lighter species. 

These results are reasonable. In single-component 
models, a Bahcall-Wolf cusp only extends outward to a 
fraction of Vm] one would not expect the pre-existing MS 
core to disappear, since it extends well beyond Tm, where 
the relaxation time is considerably longer than its value 
at Tni. Furthermore, heating by the heavier BHs should 
cause the mean density of the lighter species to decrease 



with time, as observed. 

However, at first blush, the results shown in Fig- 
ures fTSl and [Ml seem to con t radict the results described 
by iPreto fc Amaro-Seoand (|201C1( ). who used Fokker- 
Planck and A''-body integrations to follow the evolution 
of models with two mass species and a MBH. Those 
authors stated that "mass segregation... speeds up cusp 
growth [in the lighter component] by factors ranging 
from 4 to 10 in comparison w ith the single-mass case." 
IPreto fc Amaro-Seoand (J201CII ) concluded that relaxation 
to a mass segregated, coUisional steady state takes place 
in a time much less than the relaxation time at the influ- 
ence radi us, hence that coUisionally- relaxed models like 
that of .Hopman fc Alexander! ()2006D should be a good 
description of nuclei like that of the Milky Way (in spite 
of the fact that the Milky Way is not observed to contain 
a Bahcall-Wolf cusp in the stars) . 

The initial c onditi ons adopted by 

IPreto fc Amaro-Seoand (|201Q ') were rather different 
than in our post-merger models: our models have 
bona-fide cores, while their initial models had density 
cusps, p ^ r~'^ with 7 — (1/2,1). Nevertheless, the 
"acceleration" that they describe might be expected to 
occur also in our models. 

To understand the nature of this apparent discrepancy. 
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Fig. 14. — Density profiles for MS stars and BHs at different times (0.2,0.5,1,2,3) T,. during the post-merger phase. Thick sohd lines 
indicate profiles at the beginning of the post-merger phase. Dotted lines indicate r^n- 



we carried out a number of separate iV-body experi- 
ments, as well as Fokker-Planck integrations. The latter 
are presented in the next sub-section. 

We first checked whether our A^-body code would re- 
produce the rate of cusp formation observed by other 
authors in single-component models. The left panel of 
FigurefTSl shows the evolution of model CI, which had the 
same initial density profile and MBH mass as "Run 1" 
of IPreto et all ^(M \. and N = 131072. A p{r) ~ r'^/" 
density cusp forms in roughly one relaxation time at radii 
r <, 0.5rm, and the time dependence of the density is 
in excellent agreement with what was found in that pa- 
per a nd in other A^-body studies (e.g. iBaumgardt et al.l 
l2004r ). Additional experiments, varying the softening 
length, time step and particle number, convinced us of 
the robustness of this result. 

The right panels of Figure[T5] show the results of in- 
tegrating a two-component mo del, C2, with the same 
initial conditions as "Run 1" of IPreto fc Amaro-Seoani 
dlOlO), and N = 131072. Cusp growth is clearly ob- 
served only in the heavier component, similar to what 
we described above in the A^-body integrations of the 
four-component A and B. Again in this case, additional 
experiments using different integration parameters con- 
firmed the results. 



5.1. Fokker-Planck models 

We carried out integrations of the isotropic Fokker- 
Planck equation describing galaxies containing two stel- 
lar mass groups and a MBH. We used these integrations 
to address two questions. (1) What is the nature of the 
"accelerated cusp growth" that Preto & Amaro-Seoane 
observed in their Fokker-Planck integrations? (2) Is the 
evolution that we observe in our multi-component A^- 
body models consistent with the predictions of Fokker- 
Planck models? 



We begin by summarizing the evolution equations|3 
Let fi{E,t) be the phase-space number density of the 
i-th species at time t and (binding) energy E, where 
E = —v^jl -I- "0 and -0 — —'^ with $ the gravitational 
potential. Let i = 1 denote main-sequence stars, of mass 
mi, while i — 2 denotes stellar- mass BHs, of mass TO2; 
we set mijrax = 10 as in the A^-body integrations. The 
evolution equation for the i-th component is 



dU 



dF,, 



4-Mi^)^ = -^ 



dt 



F,^ 



E 

j = l,2 



dE 
Dee 



''dE 



Dsij/i 



(26a) 
(26b) 



D 



EEij 



leTT^Fm^ 



DEij — —IQTrTmim 



pE />oo 

q{E) f,{E\t)dE'+ f,(E',t)q{E')dE' 

Jo Je 

(26c) 

f,{E',t)p{E')dE' 



(26d) 



(|Merrittll20l2l) . Here, p{E) and q{E) are given by 
ri,-\E) 



p{E)^A 



r v{E, r)dr, 



"^'^U 



^-\E) 



r v'\E,r)dr 



(27a) 



(27b) 



with -if) ^(E) the inverse of the potential function, 

AttG^IhA 



v{E,r) = y^2[ij{r)-E], and F 



To 



^ These equations d if fer fr om the similar equations given by 
IPreto & Amaro-Seoand ||201(1I 1: the latter appear to be missing 
some multiplicative factors, as well as having an incorrect depen- 
dence of the diffusion coefficients on the rrij . However, we believe 
that their numerical implementation was based on the equations 
in their correct form. 
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Fig. 15. — Left; Density profiles in the single mass model CI at diff'erent times: (0, 0.25, 0.5, 1, 2) Tr(rm). The thick line indicates the 
initial profile. Right: Density profiles in the two-component model C2 at different times: (0, 0.25, 0.5, l)Tr(rm). Solid lines are from the 
Af-body integrations while dashed lines are FP models. Vertical dotted lines indicate rm in all panels. 



simplify the calculations, we made t he assumption (as 
in numerous earlie r studies e.g. IPreto et alj [20041 : 
iMerritt et alll2007al: IPreto fc Amaro-Seoanell2010D that 
the gravitational potential, 



-0(r) = 



GM, 



Air), 



(28) 



the sum of the stellar and MBH potentials, was constant 
with time, and given by its value at i = 0. This is a 
reasonable approximation at all radii: even if the stellar 
density evolves significantly at r ^ rm, the potential at 
these radii is dominated by the MBH and is nearly un- 
changing. At r ^ Tni there is no significant evolution in 
the density and the approximation is again valid. 

The coupled equations ([26| were solved by standard 
techniques, starting from initial conditions in which the 
two components had configuration-space densities 



P.W = p.(o) 



''o 



1 



r 



7 — 4 



(29) 



with the same (rg , 7 ). This is the same initia l mass distri- 
bution adopted by IPreto fc Amaro-Seoand (|2010( ) , who 
set 7 = (1/2, 1). Parameters defining the initial condi- 
tions were 7; M,/Mgai, the ratio of MBH mass to the 
total mass in components 1 and 2; and the number den- 
sity ratio 

As unit of time we adopted the relaxation time defined 
above, evaluated at the MBH influence radius rm, with 
rm defined as in the A^-body models; note that the value 
of log A becomes irrelevant when the time is expressed in 
this way. 

The right panels of Figure[Tn] show results from a 
Fokker-Planck integration of the same initial conditions 
used for the A^-body model C2 described above. The 
agreement is very good. 

Figure[TH] shows integrations of two models both with 
7 = 0.5 and M./Afgai = 0.05. The models differ in 
the fraction of heavy objects: TZ — and TZ = 10"'^. 
Th e model with TZ = 10" " ^ is th e same model plotted 
by IPreto fc Amaro-Seoand (|2010|) in their Figure 3. We 



have plotted only the density of the lighter (MS) compo- 
nent. 

The superficial appearance of these plots depends 
strongly on the radial range plotted. If the range is 
sufficiently large, extending to radii <^ rm, a MS cusp 
with slope d log p/d log r w —3/2 catches the eye in the 
model with BHs. Thi s is th e feature emphasized by 
IPreto fc Amaro-Seoand (j2010t ). Viewed more closely, in 
the region 0.02 ^ r/r^ ^ 1, the plots tell a different 
story: the model including BHs exhibits a lower density 
in MS stars at all times. Excepting at very small radii, 
addition of the BHs has the expected effect of reducing 
t he density of the MS sta rs. 

iPreto fc Amaro-Seoand (|2010t ) attributed the "accel- 
erated cusp growth" in the lighter component to "mass 
segregation," without stating explicitly the connection 
between the two phenomena. In fact, as we now argue, 
the mechanism driving the evolution of the lighter com- 
ponent at small radii is not mass segregation; it is scat- 
tering of the light component by the heavy comp onent 
(jMerritt et al.ll2007al: [Alexander fc Hopmanll2009l ). 

Consider a two-component system; as before, species 
no. 1 is the light (MS) component and species no. 2 is the 
heavy (BH) component. The four diffusion coefficients 
that appear in the Fokker-Planck equation for the light 
(MS) component scale with rrii and fi as 



DeEii 


^mjfi z 


- mipi, 


De,, 


~ m\fi ~ mipi 


DeEi2 


-TO2/2 : 


^ m2P2, 


De,2 


~ mim2/2 =i mi/92 



If 7712/02 > TTLxpi, Deei2 ^ Deeh', m othcr words, self- 
scattering is negligible compared with scattering off of 
BHs. The first-order coefficients are smaller than Deei2 
by factors of {rnipi)/{m2P2) and mi/m2 and can also be 
ignored. The evolution equation for the light component 
becomes in this limit 



5/] 



MS 



1 



d 



dt An^p dE 



D 



EE- 



9/1 



MS 



dE 



(31) 



with Dee given by equation (|26cp after setting rrij = 
tobh- The steady-state solution is obtained by setting 
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Fig. 16. — Density profiles of the lighter (MS) component in Fokker-Planck integrations of two, two-component models; the heavier 
component, the BHs, are assumed to have masses ten times the mass of a MS star. The initial conditions differ only in terms of the fraction 
of heavy particles; "R, = A^bh/^ms = (top) and 10~^ (bottom). Both components have 7 = 0.5 initially, a "core," and A/./Mg^i = 0.05. 
The model with TZ = 10~^ is the same model plotted in Figure 3 of lPreto & Amaro-Seoan3 l|201CI ) , as an illustration of "accelerated cusp 
growth;" times shown are also the same as in that figure, i.e., t = (0,0.05,0.1,0.2,0.25) in units of the relaxation time at the infiuence 
radius. Panels on the left show the MS density profile at low spatial resolution, while panels on the right focus in on the region nearer rm', 
the dotted lines on the left delineate the region p lotted on the right and the das hed lines have logarithmic slopes of —7/4 (right) and —3/2 
(left). The accelerated cusp growth described bv lPreto fc Amaro-Seoanel 120101) is seen to be present only at small radii, r < 0.05rm; it is 
due to scattering of the MS stars by the BHs, which causes the initial "hole" in phase space to rapidly fill in. At radii r > 0.05rni, adding 
the BHs has the opposite effect, resulting in a lower density of the MS component at all times. 



the term in parentheses (the flux) to zero, yielding 
dfus 



dE 



0, /] 



MS 



const 



(32) 



independent of /bh- A constant /ms corresponds, in a 
point-mass potential, to a density 

PMS cx r-3/2^ (33) 

which describes the MS cusp that appears, at early times 
and at small radii, in the Fokker-Planck models. 

The manner in which this steady state is reached will 
depend on the initial conditions. In the models consid- 
ered here, the initial MS density is pms oc r^^'^, which 
corresponds to an f{E) that tends to zero near the MBH 
- a "hole" in phase space at low energies. Scattering of 
the MS component by the BHs rapidly "fills in" this hole 
as it drives / toward a constant value. The result is a 
sharp increase in the MS density at early times at the 
smallest radii. 

The condition that the evolution of the light (MS) 
component be dominated by scattering off the heavy 
(BH) component - as opposed to self-interactions, 
which tend to build a steeper, Bahcall-Wolf cusp - is 



Pbh ^ (toms/™bh)pms ~ O.lpMS- In Figure[T71 the first 
time at which this condition is satisfied is marked by open 
circles. That figure shows evolution of the local slope, 
\d\ogpMs/d\ogr\, at two radii, O.OSrm and 0.002rni, for 
several different values of TZ, as computed via the Fokker- 
Planck equation. At the larger radius, the BHs remain a 
small fraction of the total in most of these models, and 
the evolution of the MS component is not strongly af- 
fected. But because the Bahcall-Wolf cusp builds "from 
the outside in," its initial growth is hardly reflected at 
much smaller radii. Here, as the lower panel shows, the 
dominant effect is scattering by BHs, and the effect of 
the scattering on the MS density profile is strongly de- 
pendent on (roughly proportional to) TZ. 
We summarize our findings in this section as follows. 

1 One- and two-component Fokker-Planck models re- 
produce well the behavior seen in iV-body integra- 
tions with the same initial parameters. 

2 In two-component (MS -f BH) models, addition 
of the heavy (BH) component results in a slightly 
lower MS density at radii r ^ O.OSrm, due to heat- 
ing of the stars by the BHs. 
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Fig. 17. — Evolution of the local slope, 7 = —dlogp/dlogr, 
of the MS density profile, in two- com ponent (MS + BH) Fokker- 
Planck models like those in Figure llGl The curves differ in terms of 
the heavy-particle number fraction, logjg 7?. = logj^g(A'BH/^Ms) = 
—7, —6, —5, —4, —3, —2; increasing line width corresponds to in- 
creasing TZ. Upper panel shows slopes at r = O.OSrm, bottom 
panel at r = 0.002rjn. Open circles indicate the times at which 
Pbj] = O.lpMS at the respective radii. Adding a heavy (BH) com- 
ponent has relatively little effect on the growth of a Bahcall-Wolf 
cusp in the MS component (top panel), however it greatly affects 
the form of the density profile at r < O.Olrm (bottom panel), for the 
reasons discussed in the text. The latter phen omenon is the "accel- 
erated cusp growth" described bv lPreto &: A maro-Scoane (2010). 

3 The same heating more promptly modifies the MS 
density profile at small radii, r ^ O.Obr^, at least 
in models that start from a flat core. The rate of 
growth of this "mini-cusp" is strongly dependent 
on the BH fraction. 

The fact that the effects of scattering of MS stars by 
BHs is restricted to small radii, r <^ O.OSrm, is consistent 
with the fact that we do not observe this phenomenon in 
the A^-body simulations: for instance, the density pro- 
files of FigurefHl extend down to only ^ 0.05r^. Only 
a handful of particles were ever present at smaller radii. 
At the radii resolvable by the N- body simulations, the 
Fokker-Planck models predict that the BHs should retard 
the growth of the Bahcall-Wolf cusp in the MS compo- 
nent (Figure[T6|. not accelerate it, and this is what we 
see in the A^-body simulations. 

5.2. Time dependence of the number of BHs near the 

MBH 

We now return to a discussion of the post-merger A^- 
body integrations. Figure[T3l showed the evolution of the 



100 




Fig. 18. — Cumulative radial distribution of stellar-mass BHs 
(solid lines) in Model Al. Different curves refer to diff erent times: 
(0, 0.2, 0.5, 1, 2, 3) Tr from right to left, as in FigurcfTil The dotted 
vertical line indicates rm. Dashed lines show fits to AfBH(< r) in 
the radial range [ri, r2] such that A^bh(»'i) = 5 and AfBH(^2) = 25. 




Fig. 19. — Evolution of the slope, a = — d log Nsn/d. log r, de- 
rived fro m re gression fits to the post-merger A^-body data, as shown 
in Figure lTSl for Model Al. The corresponding density-profile slope 
is 7 = 3 — o. 

mass enclosed within rm for each of the four species in 
each of the A^-body integrations. In this section, we look 
more closely at how the number of BH particles evolves 
with time on smaller scales. 

Figure[TH] plots the cumulative radial distribution of 
BHs in Model Al at different times. Time zero in this 
plot corresponds to the moment that the two, MBH par- 
ticles were combined into one. The smallest radii at 
which there are any BH particles in the A'^-body mod- 
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Fig. 20. — Number of stellar-mass BHs within 0.01 pc versus time in our models, scaled to the Milky Way nucleus, assuming (a) a density 
profile of constant logarithmic slope at small radii, and (b) a constant density inside O.lrm. 



els decreases from ^ 0.1 rm at early times to ^ O.Olrm at 
late times. In order to estimate BH numbers at smaller 
radii, we carried out regression fits of logiVBH to logr. 
The fits were performed in the radial interval [ri, r2] 
such that A^bh(< n) = 5 and A'^bh(< ''2) = 25. The 
resulting slopes, a = \d\ogNBR/d\ogr\, are plotted in 
Figure[T9] for all models as a function of time in units of 
the initial relaxation time. We find slopes 1 ^ a ^ 3 for 
NBnir), which imply slopes in the mass density profile 
in the range ^ — d log p/d log r <, 2. 

Figure[2n] shows the inferred number of BHs at r < 
0.01 pc as a function of time. In making this plot, we 
did a rough scaling of our models to the nucleus of the 
Milky Way, assuming an influence radius r^ of 3 pc. The 
factor: 



4^^0^«x!^. 3.3x103 



10 Mf; 



M. 



was used to convert the number of BH particles in the 
simulations to the actual number; the first of these fac- 
tors contains masses in physical units, the second in the 
units of the A^-body code. Since the number of BH parti- 
cles at these radii is small, we extrapolated inward under 
two assumptions: a density profile of constant power-law 
index, and a constant density inside r = O.lrm- In the 
former case, we used the slopes derived from the fits to 
-^bh(?')- Since the fitted slopes are typically large, the 
former assumption results in much larger, inferred num- 
bers of BHs. 

It is tempting to compare the numbers so ob- 
tained with estimates of Nbh obtained from steady- 
state Fokker-Planck models of the Milky Way nucleus 
(jHopman fc Alexandeii [20061: iFreitag et al.ll2006D . Here 
we note one ambiguity associated with such comparisons. 
In the Milky Way, the two infiuence radii rh and Tm are 



similar. The first is 



''h 



GM, 



3.5 pc 



M. 



4 X 106 Mq 



70 km s" 



(34) 

The second is somewhat less certain, but dynamical es- 
timates of the ma ss dist r ibutio n in the in ner few par- 
sees (e.g. Schodel et all (|2009l ): lOh et all (pOOa i) give 
Tin w 2 pc, consistent within a factor of two with r^. The 
near-equality of rh and r^ in the Milky Way is due to 
the fact that the density profile is similar to that of the 
singular isothermal sphere, p ~ r~^; the radius of the 
Milky Way's core is ~ 0.5 pc, substantially smaller than 
both rh and r^. In our post- merger iV-body models, on 
the other hand, core radii and rm are both substantially 
larger than rh. This fact precludes a unique scaling of 
our models to the Milky Way - at least in the Galaxy's 
current state. (T he Milky Way 's core may have been 
larger in the past (: Merrittll2010f ).) 

With this caveat in mind, we assume that r^ deter- 
mines the scaling of our models to the Milky Way. Fig- 
ure[20] is based on this scaling. In their c ollisionally- 
relaxed models, iHopman fc Alexandeii ()2006[ ) found 



NBuir < 0.01 pc) = 150. 



(35) 



Those authors assumed a mass function with the same 
four species as in our models, and with the same rela- 
tive numbers as in our models A2 and B2 (Table[T]). In 
Figure[20l the inferred number of BHs inside 0.01 pc is 
very uncertain at the relevant times, i.e. t w 10"'^° yr, 
fluctuating between zero, and maximum values of ^ 1 
(Model A2) and ~ 10 (Model B2). These upper hmits 
are fa ctors of ~ 10^ and ^ 10 ^ , resp ectively, smaller than 
in the lHopman fc Alexandeii (|2006l ) models. There is an 
independent way to reach a similar conclusion. After 
several relaxation times, Nbu{< 0.01 pc) is ~ 10 (Model 
A2) and ~ 100 (Model B2). These numbers, presumably 
representing the mass distribution in a near steady-state. 
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are '^10 times larger than their values at t = lO^'^ yr. 

It is interesting that the BH distribution takes so long 
to reach this (nearly) steady state - a time of at least 
twice the relaxation time as defined by the dominant 
(MS) population (Figure fTH) . This may be due in part 
to the fact that the MS distribution is also continuously 
evolving. Another reason is the persistence of a core in 
the dominant component. Chandrasekhar's (1943) dy- 
namical friction coefficient , in its most widely-used form 
(jBinnev fc Tremaind 119871 ) . predicts that the frictional 
force near a MBH drops essentially to zero if p{r) in- 
creases more slowly than r~^/^ toward the center. The 
(isotropic) phase-space density f{E) corresponding to 
that density profile has no stars moving more slowly than 
the local circular speed, and Chandrasekhar's formula 
identifies the frictional force exclusively with field stars 
moving more slowly than the massive body. This prop- 
erty of the dynamical friction force was automatically in- 
corporated into the Fokker-Planck calculations presented 
above, since that equation is based on Chandrase khar's 
coeffi cients in their standard forms (Roscnbluth et aLl 
I1957D . In reality, some part of the frictional force in the 
A'^-body simulations will come from stars moving faster 
than the test star, and including the contribution from 
these stars kee ps the frictional force fro m falling iden- 
tically to zero (jAntonini fc Merrittll201lT ). Nevertheless, 
one expects dynamical friction to act more slowly in these 
models than expected based on the application of stan- 
dard formulae that neglect the special properties of /. 

6. SHAPE AND KINEMATICS OF THE STELLAR 
SPHEROID 

We determined the shape of the merger remnants 
during the binary hardening phase by computing the 
axis ratios at different distances fro m the c enter . We 
followed the p r ocedu re described in iKatzl ()199lD and 
lAntonini et al.l (|2009[ ). We selected all particles within 
a sphere of radius d centered on the binary center of 
mass. The axis ratios were then determined from the 
eigenvalues la of the inertia tensor / as 



(36) 
where /max = rnaxjli 1,122,133}. New axis ratios were 
computed considering only particles enclosed in the el- 
lipsoidal volume having the previously determined ratios, 
i.e. all particles satisfying the condition qi < d, where 



1i ^ [^ 



(37) 



The last step was iterated until the axis ratios converged. 
If we define the axes a, 6, c such that a > 6 > c, we 
find that c/a and b/a correspond, respectively, to the 
minimum and intermediate values of ^,ri,9. From the 
axis ratios it is possible to define a triaxiality parameter 



T = 



(38) 



< T < 1, such that T = for an oblate spheroid, 
T = 1 for a prolate spheroid and T = 0.5 corresponds 
to maximum triaxiality. The ellipticity e = ^1 — (? j o? 
measures the degree of flattening of the system. The axis 
ratios, triaxiality parameter and ellipticity for models A 



and B are shown in FigureUT] as a function of distance 
from the binary center of mass, at the time when the 
binary becomes hard. A moderate triaxiality is present 
in both models at distances of the order of 1 — 2 rm 
from the binary center of mass. At larger distances, the 
models appear axisymmetric with a small flattening in 
the direction perpendicular to the binary orbital plane. 
These features persist throughout the hardening phase. 
The flattening, due to rotation of the merger remnant, 
is also visible in the isophotes shown in Figure[22l which 
are computed at the beginning of the post-merger phase. 

Rotation is introduced by the merger process, as il- 
lustrated by the velocity map in Figure!^ for Model A. 
The figure shows the direction and magnitude of the ve- 
locities in the binary's orbital plane. By the time the 
binary reaches the hard binary separation, a well defined 
rotation pattern has been established. 

The departures from spherical symmetry that we 
observe in the merger models are probably responsi- 
ble for the efficient hardening of the m assive binary 
(|Merritt fc Poon|[200l fBerczik et al.][2006[) . 

7. SUMMARY AND DISCUSSION 

7.1. (Re-)growth of Bahcall-Wolf cusps in 
multi- component systems 

In galaxies containing a single stellar population, a 
p oc r~^/"^ [Bahcall fc Wolj (J1976D cusp is expected to ap- 
pear at radii r < rsw ~ 0.2rin around the MBH. While 
the growth time of the cusp is dependent on the initial 
conditions, simulations starting from a shallow cusp in- 
side Tm show that the stellar density will have reached 
an approximate steady state after roughly one relaxation 
time at r^- We presented an example of such evolution 
in FigurclTS] 

In nuclei with two mass groups, e.g. solar-mass stars 
(MS) and 10 M0 black holes (BHs), evolution toward a 
steady state near the MBH depends on the relative num- 
bers and masses in the two groups, as well as on the ini- 
tial conditions. The results of our A''-body and Fokker- 
Planck integrations of two-component models were pre- 
sented in Section 15.11 Figure[T7] showed that addition 
of the BHs reduces slightly the rate of formation of a 
Bahcall-Wolf cusp in the MS (light) component, and also 
affects the final value of the density-profile slope, which 
varies from pMS ~ r~^'^ when pbh/pms is small, to 
-- r~^/^ as pbh/pms is increased (B ahcah fc WoHlll977[ ). 
In addition, if the initial MS distribution is very flat near 
the MBH, scattering by the heavier BHs can dominate 
the evolution of the MS component at early times for 
r < 0.05rni, converting pMS ~ r~^/^ to pMS ~ r~^^^ at 
these small radii, even before the Bahcall-Wolf cusp has 
fully formed at larger radii. 

The full galaxy merger simulations presented here al- 
lowed us, for the first time, to evaluate the evolution of 
multi-component nuclei starting from initial conditions 
that were motivated by a well-defined physical model. 
Our pre-merger galaxies contained mass-segregated nu- 
clei with four mass components, representing an evolved 
stellar population. These initial distributions were mod- 
ified both by the galaxy merger, and by the formation 
of a MBH binary, which created a large core in each of 
the components (Figures[7l [TT|) . The core radius of the 
heaviest (BH) component was somewhat smaller than the 
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Fig. 22. — Projected density contours for Model Bl at the be- 
ginning of the post-merger phase. Left: xy plane. Middle: xz 
plane. Right: zy plane. The position of the BH is indicated by 
the filled circle. The merger remnant is flattened in the direction 
perpendicular to the plane of the merger. 



Fig. 2f . — Axis ratios, triaxiality parameter and ellipticity for models A (left) and B (right) as a function of distance from the center of 
the binary, at the time when a r~j af^. The vertical lines indicate Tm. 

<^ Tr{rm), somewhat greater than in earher models wit h 
more ideahzed initial conditions (e.g. lFreitag at al.l2 006'). 

Our pre-merger galaxies (models Al, Bl) contained 
larger numbers of remnants than would be expected 
based on a standard IMF; this was done in order to bet- 
ter resolve the distributions of those components near 
the MBH. We tested the dependence of the posi-merger 
evolution on the assumed mass function by carrying out 
a second set of integrations in which we decreased the 
relative numbers of remnants (NS, WD, BH) by factors 
of a few, to values more consistent with standard IMFs 
(models A2, B2). Evolution of the dominant, MS com- 
ponent in the latter models differed only modestly from 
its evolution in the models with larger remnant fractions; 
the main difference was a lower rate of core expansion re- 
flecting a lower rate of heating by the BHs (FigurcIT^ . 
The rate of growth of the Bahcall-Wolf cusp in the MS 
component was essentially unchanged (FigurefM]). 

Our results on the regeneration of Bahcall-Wolf cusps 
following dissipationless mergers are consistent with 
tho se obtained in sim ulations of single-component galax- 
ies (jMerritt fc Szellli2006. ). We can summarize these re- 
sults by stating that regrowth of a cusp in the dominant 
stellar component requires a time comparable with, or 
somewhat longer than, the relaxation time of that com- 
ponent measured at the MBH influence radius. The new 
simulations presented here suggest that time scales for 
cusp regrowth are only weakly dependent on the number 
of heavy remnants (BHs) , at least if the fraction of mass 
in the heavier population does not exceed a few percent 
of the total. 

The Milky Way nucleus is near enough that a Bahcall- 
Wolf cusp could be resolved if present, and the dom- 
inant stellar population is believed to be old. Since 
?'h ~ ''m ~ 2 — 3pc in the Milky Way, the expected, outer 
radius of the Bahcall-Wolf cusp is rew ~ 0.5 pc w 10". 
As is well known, number counts of the dominant, old 
stellar population show no evidence of a rise in density 



cores in the three lighter components, a relic of the earlier 
mass segregation, and of the incomplete cusp destruc- 
tion process during the binary MBH phase (Figure fTTj) . 
Evolution during the merger phase set the "initial con- 
ditions" of the nucleus at the time when the two MBHs 
were combined into one. Unlike the rather ad hoc initial 
conditions used in man y earlier stu dies of nu clear evo- 
lution (e.g. iFreitag et a l. 2006; Pre to fc Amar o-Seoane 
l201Clf ). our post-merger models have cores that should 
be reasonable representations of the cores in real galax- 
ies that formed via dissipationless mergers, with mass 
densities that decline toward the center and anisotropic 
kinematics that reflect the action of the massive binary 
on the stellar orbits (Figure E5|) . 

By continuing the evolution of these multi-component 
models for a time greater than Tr(rm) after coalescence 
of the two MBHs, we found that the cores characterizing 
the distribution of the dominant, MS component per- 
sisted; in fact the MS density at radii ^ r,„ decreased 
gradually with time (Figure 1131) - a predictable conse- 
quence of the initial extent of the cores, several times 
Tm, implying Tr{rc) > Ti.(r,„), and of continued "heat- 
ing" by the heavier BHs. Nevertheless, a Bahcall-Wolf 
cusp gradually reformed in the lighter components at 
radii r <^ 0.2rni <?C Tc', growth times were found to be 
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Fig. 23. — Velocity vectors in the orbital plane (x ~ y plane) of the binary at different times during the merger of Model A (i = 
50, 125, 250, 375, 425, 625, from the top left to the bottom right). The origin corresponds to the initial position of the binary center of mass. 
The velocity vectors represent the in-plane velocities of all stars in the chosen area. The circles represent the locations of the massive black 
holes. The final frame is shown at higher resolution. 



at this radius; instead the number counts are flat, or 
even falling, from ^ 10" into at least 1" pro.iected ra- 
dius (jBuchholz et al.ll2009t iDo et ani2009t iBartko et alJ 
[2OIOI) . 

Assuming solar- mass stars, the relaxat ion time at th e 
influence radius of Sgr A* is 20-30 Gyr fMe rrittilmOl ). 
while the mean stellar age in the nuclear star clus - 
ter is estimated to be ■^ 5 Gyr (jFiger et al.l 120041 ) . 
The time available for formation of a cusp is therefore 
(5/25)T^(r„i) « 0.2T^(rm). Our simulations (e.g. Fig- 
ure lTil) suggest that a Bahcall-Wolf cusp in the dominant 
component is unlikely to have formed in so short a time, 
and this is consistent with the lack of a Bahcall-Wolf cusp 
at the Galactic center. 

The core in the Milky Way has a radius of ~ 0.5 pc, 
somewhat smaller than rh or rm, while the cores formed 
in our merger models are somewhat larger than rm (Ta- 
ble 4). It has been argued that the Milky Way core is 
small enough that gravitational encounters would cause 
it to shrink appreciably in 10 Gyr, as the ste llar dis- 
tribu tion evolves toward a Bahcall-Wolf cusp (jMerrittl 
|2010( ) . The cores in our A^-body models are so large that 
they do not evolve appreciably after the binary MBH has 
been replaced by a single MBH; the Bahcall-Wolf cusp 
forms at radii smaller than r^^ leaving the core structure 
essentially unchanged. Without necessarily advocating 
a merger model for the origin of the Milky Way core 
(it is unclear whether our galaxy even contains a bulge; 
[Martinez- Valpuesta fc GerhardI (|2011[ )). we note that the 



sizes of cores formed by binary MBHs scale with the bi- 
nary mass ratio. Presumably, we could have produced 
cores more similar in size to the Milky Way's if we had 
adjusted this ratio. 

It has been suggested (jLockmann et al.l I2010D that 
heating by BHs could be responsible for the lack of 
a Bahcall-Wolf cusp at the Galactic center. We do 
not find support for this hypothesis, either in our four- 
component A^-body models, nor in our two-component 
Fokker-Planck models. The latter showed that even a 
quite large BH population still allowed a Bahcall-Wolf 
cusp to form in the MS stars on a time scale of ~ Tr(r„i); 
the main effect of the BHs is to decrease the asymptotic 
slope of the cusp from ~ r"^/"* to ^ r""^/^. 

7.2. The distribution of massive remnants in galaxy 

nuclei 

A dense cluster of stellar-mass BHs has been invoked 
as a potential solution to a number of problems of coUi- 
sional dynamics at the Galactic center. Examples in- 
clude randomization of the orbits of yo ung stars via 
gravitational scattering (e.g. IPerets et al.i r2009) , produc- 
tion of hype r- velocity stars through encounters with BHs 
(e.g. lO'L ea rv fc Locb 2008), and warping of young stellar 
disks (jKocsis fc Tremaindl2011l) . These treatments typ- 
ically assume a coUisionally-relaxed state for the Galac- 
tic center. In the relaxed models, the mass in BHs in- 
side 0.1 pc is ~ IO^'Mq, i.e. A^bh (?- < 0-1 pc) w 10=^ 
and NBii{r < 0.01 pc) « 10^ (Hopman fc Alexandeil 
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120061 iFreitag et all [20061) . assuming "standard" IMFs. 
A high density of stehar BHs at the centers of galaxies 
hke the Milky Way is also commonly assumed in discus- 
sions of the EMRI (extreme - mass-ratio inspiral) problem 
(jAmaro-Seoane et an i2007'). 

Models like these are called into question by the lack of 
a Bahcall-Wolf cus p in the late-type stars at the center 
of the Milky Way d&iic hholzet aLl l2"0Q9l: iDo et al.ll2009l: 
iBartko et al.ll2010() . If the Galactic center is not collision- 
ally relaxed, as these observations seem to suggest, then 
computing the distribution of the heavy remn ants be- 
come s a more difficult, time-dependent problem (jMerrittl 
l2010f) , and knowledge of the initial conditions is essential. 

Dissipationless mergers imply "initial conditions" char- 
acterized by a core in the dominant stellar component. 
The cores formed in our (major) merger simulations are 
large enough that they do not evolve appreciably (i.e. 
shrink) even after ~ 3 relaxation times at the MBH influ- 
ence radius. As a result, evolution of the BH distribution 
takes place against a stellar background with a very dif- 
ferent density profile than in the steady-state models. In 
a core around a MBH, dynamical friction is much weaker 
than one would estimate by plugging the local density 
into standard formulae for orbital decay, due to the ab- 
sence of stars moving more slowly than the local circular 
velocity (|Antonini fc Merrittir2011[ ). Scaling our iV-body 



models to the Milky Way, we predicted numbers of BHs 
inside r^ that were substantially smaller than in the col- 
lisionally relaxed models; in the inne r 10~^ pc - the radi i 
most relevant to the EMRI problem ([Merritt et al.ll2010t) 
- the number of BHs was, at most, 10-100 times smaller 
than predicted by these models, even after 10 Gyr (Fig- 
ure l^ . 

It is unclear how relevant merger models are to the cen- 
ter of the Milky Way. If the core observed at the Galactic 
center had some other origin, the distribution of stellar 
BHs might have little connection with the distribution of 
the giant stars. However, if cores of size Tc ~ rh are com- 
mon features of galactic nuclei, and if at some early time 
both the BHs and the stars had a common core radius, 
our models imply that the distribution of BHs should be 
considered very uncertain, even in galaxies for which nu- 
clear half-mass relaxation times are as short as the age 
of the universe. 
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